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ABSTRACT 



The Marine Gas Turbine control systems in present use in the U. S. 
Navy are of such significant technological age that new design 
techniques could lead to more optimum performance and increased 
plant efficiency. To this end, a new real time Marine Gas Turbine 
simulation method is needed for advanced controller design and 
implementation. 

A modeling method is shown which utilizes real time sequential 
linearizations to approximate the true nonlinear response of the NPS 
Boeing 502-6A test facility. A validation of this simulation approach is 
presented. The method has immediate application to advanced 
controller design, especially to the design of modem regulators 
(Linear Quadratic), model reference controllers, and real time 
diagnostics. 
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I. INTRODUCTION 



Conventional control theory is limited to time invariant single-input- 
single-output systems and is usually based on a frequency domain 
design approach. Modem control theory, on the other hand, can be 
well applied to multiple-input-multiple-output, non-linear gas turbine 
systems which are now in use and it relies on a much more applicable 

time domain approach. 1 The feasibility of modern control theory 

such as Linear Quadratic Regulator (LQR) theory, in the application of 

marine gas turbines for propulsion has been demonstrated by past 

work at the Naval Postgraduate School in Monterey, Califomia.2-4 For 

this application to be successful a simple, accurate, and robust real 

time marine gas turbine simulation must be developed. 

This thesis introduces a modeling approach to be applied in 
developing a real time simulation of engine response to meet the 
requirements of simplicity, accuracy, and robustness. It is simple in 
that it is only first order dependent, accurate in that it provides 
acceptable (+/- 10%) response for use in application of advanced 
engine control design, and robust in that it works well over the entire 
operational range of the gas turbine system. 
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n. PAST WORK EN GAS TURBINE MODELING 



Past work in Gas Turbine Modeling is discussed in terms of two 
areas: first, a brief discussion on gas turbine modeling in general is 
presented. This is followed by a discussion on the previous work 
completed at the Naval Postgraduate School in Monterey, California. 

Since most marine gas turbines have been derived from aviation 
sources, most of the early gas turbine modeling work was 
accomplished in the aviation arena. In the 1970’s the Naval Gas 
Turbine Ship Propulsion Dynamics and Control Systems Research and 
Development Program was started to address the issues of designing 

Marine Gas Turbines for future use in the Navy. 5 This included the 

dynamic modeling of the response of a marine gas turbine. Prior to 

this program, most of the work had been steady state modeling. 

However, the Navy realized that dynamic modeling was needed for 

advanced controller design. The program was successful in developing 

a gas turbine machinery dynamics and control system data base. From 

this work a computer-based propulsion control testbed was developed 

for use in gas turbine controller development. Cut and try linear 

control design then ensued, which resulted in a gain adjusted 
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controller design in use in today's fleet. This approach resulted in a 
linear Proportional-Integral theory being applied to a nonlinear 
system. The results were deemed adequate for the Navy's uses. 

A new approach to this modeling is to develop a method to 
accurately model the propulsion plant nonlinear responses with a real 
time simulator. This would help eliminate or alleviate problems 
associated with the varying propeller loads felt in the gas turbine by 
predicting the responses before they actually occur, thus allowing 
more acurate control. More accurate control would mean better 
engine response to these loads and would, in turn, mean less wear, 
thus leading to longer engine life as well as better engine fuel 
efficiency. This is the direction taken at the Naval Postgraduate 
School, (NPS), in the gas turbine modeling problem. 

The previous work started with Phillip Johnson in 1985 in which 
he conducted the development and implementation of a load control 
system ( e.g. a water brake dynamometer) to emulate the shaft and 
propeller of a marine gas turbine.^ The next step was taken by James 
Roger also in 1985. He used a simple multiport diagram and the 
Continous System Modeling Program (CSMP III) to incorporate an 
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improved model for the dynamometer unloading valve into the existing 

system dynamic model. 7 These two works formed the base for the 
next step performed by Vincent Herda in 1986. Herda developed a 
better multiport model and used it to develop a linear state space 

model and a nonlinear dynamic propulsion model. 2 Robert Miller 
followed in the same year using Herda’s models and examined the 
feasibility of using Linear Quadratic Regulator (LQR) control design 
theory for a small gas turbine, specifically the Boeing 502-6A test bed 

installation at NPS.3 Finally Vincent Stammetti in 1988 developed the 
first real time dynamic simulator for the gas turbine and compared an 
LQR controller to a classical Proportional Integral (PI) regulator based 

on his results with the simulator. 4 To summarize, Herda validated 
cause and effect through the use of his new multiport model, shown in 
figure 1 , Miller validated the steady-state model method and 
Stammetti followed with the first linear model structure. It is 
important to note that these models were only accurate for one run 
of simulation and were not applicable over the whole range ol 
operation of the gas turbine. Their main point however, was to 
validate and show the feasibility of developing a better model for use in 
modern controller design to update those in current use today. To 
this end they were successful. This paper will continue this work and 
attempt to develop a simple, robust, and accurate real time marine gas 
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Figure l. Herda's Multiport Diagram. 



turbine simulator to be used for advanced controller design. 
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HI. MODEL STRUCTURE AND COARSE TUNING 



A. Modeling 

The gas turbine plant to be modeled was the 175 horsepower 
Boeing 502-6A which was assembled in a test bed configuration at the 
Naval Postgraduate School in Monterey, California. It was coupled to a 
Clayton 17-300 water brake dynamometer which was used to simulate 
the shaft and propeller loading on a marine gas turbine. The gas 
turbine /water brake assembly is shown in figure 2. The gas turbine 
itself was composed of two aerodynamically coupled sections, a gas 
generating section and a power output section. The gas generator had 
three major components. These were the compressor, the burner, 
and the high pressure turbine shown in figure 3. Connected by 
gearing off the gas generator shaft was the accessory section 
consisting of the fuel pump, lube oil pump, governor, tachometer 
generator, and the starter. The compressor was a single stage 
centrifugal compressor and was coupled to the high pressure turbine 
aerodymanically via two through-flow type, cross-connected 
combustion chambers. The high pressure turbine was a single stage 
axial flow high pressure turbine. It was aerodynamically connected to 
the low pressure or free power turbine which was also an axial flow 
turbine. The free power turbine output shaft was mechanically 
coupled to the water brake dynamometer which absorbs the energy of 




175 Hp Boeing 502-6A gas turbine 

Single stage centrifugal compressor 
Two single stage axial turbines 



Clayton 17-300 water brake dynamometer 



Assorted sensors 
Pressure 
Temperature 
Output torque 
Rotational speeds 
Fuel flow rate 
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Figure 3. Gas Turbine Components. 



t±ie free power turbine in simulating the action of a drive shaft and 
propeller. 



The modeling began with Herda’s cause and effect multiport 
diagram shown in figure 1. The importance of this diagram to the 
method can hardly be overstated. Each of the components, the 
compressor, high pressure and free power turbines, and the gas 
generator and dynamometer shafts, were modeled by describing their 
outputs in terms of their inputs in equation form. For example, the 
compressor inputs are gas generator speed (Ng) and the high pressure 
turbine inlet pressure (P2) and the outputs are mass flow rate of air 
(Ma), high pressure turbine inlet temperature (T2), and compressor 
torque (Qc). The compressor block diagram is shown in figure 4. 

Ma 
T2 

P2 Figure 4. Compressor Block Diagram. 




It is important to note that these input/output quantities as well as 
those for the other components must be measurable or calculable. 
They were recorded and or calculated over the operating range of the 
gas turbine (shown in figure 5), and the results are tabulated in table 
1. For clarity, the operating region of the turbine is shown in figure 5. 
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Figure 5. Gas Turbine Operating Range. 
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(A) Raw data 
(B) Normalized data 

Table 1. Steady State Data Taken Over Operating Range 
(Ng 20k -32k rpm, Ns 500 - 2500 rpm) 



1 2 



At this point it should be noted that the data of table 1A was 
normalized according to each of the variable's midpoint for ease of 
understanding each variable's relative impact with respect to the 
others. Once normalized, the variables were are denoted using lower 
case letters. The resulting equation forms are: 



Ma = fi (Ng, P2), 




(1) 


T2 = f 2 (Ng, P2), 


and 


(2) 


Qc = f 3 (Ng, P2). 




(3) 



The normalization values are shown in table 2. 



Table 2. Normalization Points. 



Ng = 26000 
Qd/Qfpt = 220 
P2 = 14 
P4 = 1.9 
Ma = 2.36 



Ns = 1500 
T2 = 215 
T4 = 956 
Mf = 122.5 
Qc/Qhpt = 80 



The next step was to distinguish between steady and dynamic 
components. Once this is complete, the component modeling can 
proceed. The rest of the component equation forms are shown in 
table 3. 



Table 3. Component Equations. 



Gas Generator Shaft: Ng = f* (Qhpt, Qc) (4) 

Dynamometer Shaft: Ns = fs (Qfpt, Qd) (5) 

Free Power Turbine: P4 = f6 (Ns, T4, Maf) (6) 

Qfpt = {7 (Ns, T4, Mai) (7) 

High Pressure Turbine: E = fs (Mi) (8) 

Maf = f 9 (Ma, Mi) (9) 

T4 = fio (Ma, T2, P4, Ng, E) (10) 

P2 = fn (Ma, T2, P4, Ng, E) (1 1) 

Qhpt = fi 2 (Ma, T2, P4, Ng, E) (12) 



B. State Selections 

State variables are the smallest set of variables which are required 

to ascertain the state or condition of a dynamic system. 1 The 
knowledge of these variables at an initial point in time (to), when 

combined with the input as time proceeds (t > to), will completly 

describe the systems behavior. The state of a system at any time (t) is 
indepentent of any prior state and input which precede the initial 

time point (to) ^ . Here, the turbine and load is the dynamic system and 
the states, when properly applied, will accurately describe the 



dynamic behavior of the system. Significant prior research has 
determined that only three states are necessary to sufficiently 

describe the operating condition of the NPS gas turbine.® In reality 
however, only two states are necessary (ng, and ns) since the third 
state (e) represents a short first order lag of the input fuel (mf). This 
time lag is due to the burner action delay between the input of 
additional fuel and the thermodynamic energy output felt in the 
turbine. 

It is our goal to relate all of the dependent variables (p2, 
p4,t2,mf,...) to the three state variables ng, ns, and e to allow us to use 
the state space equation x = Ax + Bu (ref. 9) to describe the transient 
behavior of our system. Here x is the state vector and is composed of 
state variables ng, ns, and e, u is the input vector and is composed of 
mf and qd. Refering back to figure 1, note that we have redrawn 
Herda's system boundary at qd in order to have a more measurable 
input. Again, from prior research, it was known that each of the states 
is associated with a significant time lag in the system response. In 
this way, we knew that the fuel input component, the rotor shaft, and 
the output shaft must be regarded as dynamic components. We also 
knew that the remaining components could be regarded as 
instantaneous (memoryless), or steady components. The steady and 
dynamic components are shown in table 4. 



Table 4. Steady and Dynamic Components. 



Steady Components: 

compressor ma 

high pressure turbine p4 
power turbine t4 

Dymamic Components: 

rotor shaft ng 

load shaft ns 

burner lag e 



t2 qc 
qfpt 

p2 qhpt maf 



C. Steady Components 

While the steady components of the model were given above as non 
linear functions, in this section, we reduce these to linear 
expressions, where the coefficients were regressed and varried 
according to the engine operating point. 

The steady component equation form selection was accomplished 
through the use of an equation fitting program (Minitab), and plotting 
of each of the input variables versus an output variable minus the 
residual. If the residual plot displayed a parabolic shape, it was an 
indication that the input variable function in question was of a higher 
order. In the case of the equations with only two input variables a 
different method was used to fit the data. That is, one at a time, each 



of the input variables were held constant and the other was plotted 
versus the output variable. Both methods are developed further in 
Appendix C. This resulted in accurate data fit equations. 

To simplify the equation reduction methodology, the steady 



component equations were kept in the form: 

x = cl y + c2 z (13) 

For example, the compressor equations appeared as follows: 

t2 = cl ng 2 + c2 ng + c3 p2 + c4 (14) 

ma = c5 ng 2 + c6 ng + c7 p2 + c8 and (15) 

qc = c9 ng 2 + clO ng + ell p2 + cl2 (16) 



where the coefficients cl thru cl 2 were the values from the equation 

fits. It is important to note the strong dependency of the state ng. 
This will appear again in the decoupling equations developed at the 
end of this chapter. The next phase of development was the 
formulation of the dynamic component models. 

D. Dynamic Components 

The dynamic components are derived from the gas generator shaft, 
the dynamometer shaft, and the lag relationship between mf and e. 
The dynamic components were modeled with differential equations 
which, when integrated over time will determine the state of the 



system. The dynamic components took the form: 





ng = qhpt-qc 


(17) 




ns = qfpt-qd 


(18) 


and 


e = (mf-e)/T. 


(19) 



These, along with the previous equations of the steady components, 
were next linearized to form a dynamic equation set for the simulation 
model. 

E. Linearized Dynamic Equation Set 

The linearized dynamic equation set was developed from the 
perturbational or derivative values of the steady state equations 
(symbolically, the perturbational variables will be proceeded by a "d"). 
The compressor equations are used below to demonstrate the 
formulation of the dynamic equations. Recall that the t2 steady state 
equation was: 

t2 = cl ng 2 + c2 ng + c3 p2 + c4 (14) 

So, the dynamic equation took the form: 

dt2 = [ 3t2/8ng ] dng + [ 8t2/8p2 ] dp2 (20) 

Here the partial differentials will be denoted with upper case letters 
so that the resulting form of equation 20 is: 
dt2 - Ai dng + A2 dp2. 
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( 21 ) 



The rest of the steady and dynamic equations previously developed 
were also converted to dynamic form in this manner. The results 
appear below in table 5. 

Table 5. Dynamic Equation Set. 

Compressor dt2 = A1 dng + A2 dp2 (22) 

dma = Bl dng + B2 dp2 (23) 

dqc = C 1 dng + C2 dp2 (24) 

Gas Generator Shaft 

JGB dng = dqhpt-dqc (25) 

Dynamometer Shaft 

JDB dns = dqfpt-dqd (26) 

Free Power Turbine 

dp4 = Dl dns+D2 dt4+D3 dmaf (27) 

dqfpt = Fl dns+F2 dt4+F3 dmaf (28) 

High Pressure Turbine 

de = (dmf-de)/T (29) 

from continuity— dmaf = Kl dmf + K2 dma (30) 

dt4 = Gl dma + G2 dt2 + G3 dp4 + G4 dng + G5 de (31) 

dp2 = HI dma + H2 dt2 + H3 dp4 + H4 dng + H5 de (32) 

dqhpt = U dma + 12 dt2 + 13 dp4 + 14 dng + 15 de (33) 



These resulting equations from table 5 will be used to develop the 
state space relationships through a reduction methodology to get 
three equations ( one each for ng, ns, and e) in terms of ng, ns, e, mf. 



and qd. 



F. Decoupling the Dynamic Equation Set 

In the past, during the equation set reduction, a number of the table 
3. equations were used more than once. This resulted in the 
development of singularity points in the simulations. To avoid this, 
another method was developed in the present work. This was 
achieved by inspecting the system block diagram (fig. 1) and 
deciding that if two equations, one for p2 and the other for p4, 
could be expressed in terms of the states ng, and ns only, then it 
would be possible to reduce the equation set directly (without 
resubstitution) thus eliminating the singularities. The new p2 and p4 
equations were developed with the assistance of the Minitab 
regression routine and steady state data for ng and ns. This resulted 
in the conclusion that both variables had a second power effect with 
ng and a first power effect with ns. With these relationships in mind, 
the computer routine was used to fit the data in a two step process. 

First , equations for p2 and p4 were generated in terms of ng A ^, ng, 
and ns. Next, the data generated from these equations for p2 and p4 
was plotted versus ng at five different constant ns values. These plots 
were compared to the actual data to check for accuracy of the 
equations. The results are shown on the two graphs in figure 6 and 
show the robustness of these equations over the entire range of the 
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P4 (PSIC) P2 (PS1G) 

0.7 1.2 1.7 2.2 2.7 3.2 0.0 1 1.0 10.0 21.0 



P2 VS NG (NS=500-2500 RPM) 




P4 VS NG (NS=500— 2500 RPM) 




Figure 6. P2 and P4 graphical Representations. 
Note: Solid lines are equations, symbols are data. 
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data. The new dynamic equations for p2 and p4 are now in the form of: 



p2 = cl ng A2 + c2 ng + c3 ns + c4 and (34) 



p4 = c5 ng A2 + c6 ng + c7 ns + c8. 
These equations were then linearized into the form: 



(35) 



dp2 = Hi dng + H2 dns for eqn. 11 and (36) 



dp4 = Dl dng + D2 dns for eqn. 6. 



(37) 



These twelve equations were then reduced with the reduction 
methodology as described in detail in Appendices A and B. 

The results of the reduction transformed the linearized equations to 
state space equations, which were used in the simulation model. The 
resulting state space equations from Appendix B are: 



where Z1 through Z7 were derived from combinations of the A1 
through K2 coefficients of the dynamic equations. The standard state 

space equation is x = A(t) x +B(t) u.9 Here x is a matrix made up of 
dng, dns, and de, x is a state vector composed of the state variables: 
dng, dns, and de, and u is a vector composed of inputs dmf and dqd. 



dng = Z1 dng + Z2 dns + Z3 de 
JGB JGB JGB 



dns = Z4 dng + Z5 dns + Z6 de + Z7 dmf - 1 dqd 
JDB JDB JDB JDB JDB 



de 



_L_de + 1 dmf 
T T 
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The A and B matrices are: 



A13 = Z2l 
JGB 



All = Z1 
JGB 



A12 = Z2 
JGB 



A21 = Z4 
JDB 



A22 = Z5 
JDB 



A23 = Z6 
JDB 



A31=0 A32 = 0 A3 3 = -1 

T 



and 



Bll = O B12 = O 



B21 = ZZ 
JDB 
B31 = J_ 
T 



JDB 
B32 = O 



B22 = zl 



and the state space equations now appear as: 

dng = A1 1 dng + A12 dns + A13 de + Bl 1 dmf + Bl2 dmf 
dns = A2 1 dng + A22 dns + A23 de + B2 1 dmf + B22 dmf 

de = A31 dng + A3 2 dns + A3 3 de + B31 dmf + B32 dmf. 

A coarse tuning was next performed to determine the best 
numerical values of the coefficients. This was necessary since the first 
set of regressed coefficients did not provide accurate simulation. This 
was attributed to the coarse instrumentation that is present on the 
NPS turbine and the manner of fitting in the regression routine. 
Consequently, attention was focused on the following coefficients as in 
need of tuning: Al, A2, Bl, B2, Cl, C2, FI, and II - 15 . Those for Dl, 
D2, HI, H2, Kl, and K2 were found previously using the data 
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regression program. Initially chosen for a starting point was that for 
B2 (3ma/3p2). This was an obvious choice for much is known about 
the relationship between ma and p2 in compressors. To find B2, p2 
was plotted versus ma at constant ng values. From this the slope 
3ma/3p2 (B2) was obtained. The B1 expression (3ma/3ng) was found 
by regression using the plotted value of B2 in the Minitab routine. A2 
and C2 were then found using plotted values to find the slopes 3t2/3p2 
(A2) and 3qc/3p2 (C2). These were then used as inputs into the 
relationships: 

B1 = 3ma/3n g = A1 = 3t2/3ng = Cl = 3qc/3ng 

B2 3ma/3p2 A2 3t2/3p2 C2 3qc/3p2 

from which A1 and Cl were obtained. The results were then checked 
and verified through the use of the data regression program. At steady 
state qc = qhpt is required, thus the relationship for Cl (3qc/3ng) = 14 
(3qhpt/3ng). So, with 14 determined, the value of II (3qhpt/3ma) can 
be found as II = I4/B1 and 12 (3qhpt/3t2) is I4/A1. With these three, 
II, 12, and 14, as inputs, 13 (3qhpt/3p4) and 15 (3qhpt/3mf) were 
obtained through the use of the data regression routine. The 
numerical values can be found in Appendix C. 

The next step involves fine tuning the computer simulation model 
and results in the setting of the remaining coefficients of the A and B 
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matrices. 



IV. FINE TUNING THE MODEL 



In line tuning, the actual dynamic data was used in conjunction with 
a simulation approach to select final values for the model constants. 

The simulation approach used the simple linear state space 
equations for dng, dns, and de. These described a first order 
differential equation set, with time varying A and B matrices, to mimic 
the actual nonlinearity of the system response. This was accomplished 
by approximating the true nonlinear relationships with a series of 
small linear steps as depicted in figure 7. 

The dynamic simulation computer model was developed from a 
routine written by Herda and later modified by Stammetti, it is shown 
in figure 8. The major change to the previous program was in the A 
and B matrix equations. In the past, the equations were complex 
and were a combination of exponentials and second order equations. 
In the present model they were simple, first order linear functions of 
the states. Another change which occurred in the fine tuning process 
was the adjusting of the dynomometer shaft inertia. This change 
was for shaping the response and did not affect the steady state 
response of the model. 

A. Computer Simulation 

The simulation program was built using four sections. These were 
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Input: 



A(t), B(t) 



Compute: 



x = A(t) x + B(t) u 



-f 



x dt 



Figure 8. Simulation Process. 
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the dynamic data input, the initial condition, the dynamic, and the 
derivative sections. The actual dynamic data was taken on a four 
channel strip chart recorder and depicts the real turbine response of 
Ng, Ns, Qd, and Mf as a function of time. These dynamic runs were 
compiled in Appendix E and provided the source and data for the 
simulation. The initial condition section established the initial 
conditions for each particular run. The dynamic section formulated 
the display of the dynamic data and was followed by the derivative 
section which formulated the simulated response. The outputs of 
these final two sections were then plotted for comparison. These 
plotted results were used to assist in a grid search on the still doubtful 
component coefficients to determine the final values of the A and B 
matrices. The procedure was to guess values of the coefficients, then 
simulate the model response. If the model and the data agreed, then 
the search was over. 

B. The Grid Search 

A grid search technique was applied in two and three parameter 
searches to find the optimal values for the remaining coefficients. It 
was formed using a the spread sheet formatted program called Excell, 
by Microsoft software. An example of this is found in Appendix C. The 
top portion provided the inputs for the coefficients used in the 
calculation of the A and B matrices. The middle portion calculated the 
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intermediate Y variables. These were then combined with the 
previously mentioned A thru K coefficients to calculate the A and B 
matrix values, three of which were ng dependent. 

We also knew that the steady state required x = 0. This fact was 
used to generate additional equations which must be satisfied by any 
coefficient set to guide the grid search. The last section of the spread 
sheet provided for the calculation of the values of dng and dns for the 
simulation runs 3 and 9 as well as the peak dns value for run three 
(figure 9). Run three, shown in figure 9, developed into the most 
difficult run for determining the best coefficients for Ns simulation 
and is the reason for the addition of the peak dns value. The values 
for the G, F2, and F3 coefficients were the focus of the grid search as 
we had least confidence in them. These are the important coefficients 
in determining the value of dns and in turn, the ns simulation values. 
The G1 - G5 values were calculated in the same manner as described 
eariler for the 11-15 values. They represent the partial derivatives of 
the equation for t4 (eqn. 10) and were determined to have little effect 

on the A21, A22 and A33 results. Their elimination as a search 

parameter reduced the search to a two parameter one. Various values 
for F2 and F3 were tried, and the resulting A and B matrix coefficients 
calculated were inserted in the computer simulation. The results 
were plotted on the NPS mainframe computer system, using the 
easyplot routine, to check for accuracy of simulation. This process was 
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32000 

29000 

ng 

26000 

(RPM) 

23000 

20000 
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500 1000 1500 2000 2500 

Ns (RPM) 

Figure 9. Dynamic Data Runs. 



repeated until the optimal F2 and F3 values were obtained. 
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V. VALIDATION OF THE MODEL 



The validation of the model was a two step process which carried 
over from the fine tuning. After the fme tuning was complete, and 
the model resulted in the proper shape for runs 3 and 9, and with the 
steady state response within the required tolerances, the first step 
was complete. The next step was to take the model and apply it 
elsewhere in the operating range of the turbine to check for accuracy. 
If the proper response in shape and prescribed tolerance criteria (+/ 
- 10 percent in the dynamic and steady response) were still met. 
then the validation was completed. This process was only performed 
around the borders of the operating range (figure 9). If successful, 
then the model should be applicable anywhere in the turbine 
operating range. 

The first run modeled was the one for run nine which models the 
Ns transient response from 500 to 2500 rpms while holding the value 
for Ng, as load varies, constant at 32000 rpms (+/- 5%) for the 
conditions shown in run nine in Appendix E.. Figure 10 is the 
graphical representation of the modeled run and the dynamic data run 
for this transient. The criteria has been met for run 9, so the model 
was applied to the validation run one. Run one, simulates the 
transient for Ns from 500 to 2215 rpms (the turbine will not operate 
above 2215 rpm at low Ng) with the value for Ng constant at 20000 
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NS (RPM) NG (RPM) 

500 1000 1500 2000 2500 3000 20000 23000 20000 20000 32000 35000 



SIMULATION RUN 9 NG 




TIME (SEC) 



SIMULATION RUN 9 NS 




Figure 10. Simulation-Run Nine. 
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NS (RrM) 

500 1000 1500 2000 2500 3000 NG (RPM) 



VALIDATION RUN 1 NG 




TIME (SEC) 



VALIDATION RUN 1 NS 




Figure 11. Validation-Run One. 
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rpms (+/- 5%); the results are shown in figure 11. Again, the 
accuracy criteria were met. This covers the upper and lower bounds 
of the operating range. 

Next, run three was used in conjunction with run nine in the fine 
tuning portion. This run is a vertical run on the left side of the 
operating range and is an Ng transient from 20000 to 32000 rpms 
with Ns being held constant at 500 rpms to conform with the actual 
dynamic data run. This run is shown in figure 12. Note that the 
accuracy criteria are easily met in Ng but apparently not satisfied in 
Ns. However, when the accuracy criteria is applied to the normalized 
Ns, the transient is acceptable. To validate this run the model is next 
applied to the right hand side of the operating range. This is run 
seven and, again, is a Ng transient from 20000 to 32000 rpms with Ns 
slightly increasing from the lower constraint of 2215 rpm to 2500 
rpms which is attainable at the upper limits of the operating range. 
This run is presented in figure 13. A validation run was attempted 
vertically in the center of the operating range. This is run five and 
it's response was of the same result shown in figure 1 1 for run seven 
and therefore was not shown as a separate figure. 
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Figure 12. Simulation -Run Three. 
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VALIDATION RUN 7 NG 




VALIDATION RUN 7 NS 




Figure 13. Validation-Run Seven. 
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VI. CONCLUSIONS 



A modeling method has been developed for real time simulation of a 
marine gas turbine. The results of this modeling method should be 
used with the LQR controller design theory to develop a new 
controller for the Boeing test bed at NPS. This controller should then 
be tested to ascertain it's ability to be more fuel efficient than the 
present system. 

This model was simple in that it used matrix elements for A and B 
which were only linear state functions. The dependency was only on 
the Ng state. This simplicity will provide fast run times on a small 
computer for model reference control and diagnostics. 

Accurate simulation performance has been demonstrated in the 
compressor/high pressure turbine section and the resulting steady 
states have been achieved to +/- 3%. An exception to this accuracy is 
the constant Ns simulation for the vertical runs . This is only an 
indication that the final values for F2, F3, and possibly G1 are not fine 
tuned enough to provide for more accurate constant Ns response. 
Investigation was conducted into this irregularity to determine the 
cause of the steady state results being to high. The fact that the steady 
state response for a constant Ns is higher in runs five and seven was 
probed and a number of solutions were looked at. One such solution 
was to check for the possibility that the oscillating torque for runs five 
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and seven shown in Appendix E was causing the Ns value to continue 
to rise and settle out at a higher steady state. These oscillations do 
not appear to be the cause nor have any affect on the steady state 
response of Ns. This was verified by modifying the computer 
simulation routine so as to remove the oscillation. Another solution 
was to adjust the A21 values until a proper steady state value for ns was 
reached. At the time of publication this has not provided any 
satisfactory results. Another possibility is that the A21 coefficients 
may have to change the sign to produce the proper response. The last 
area to examine as a possible solution to this problem is to input a 
varying inertia for the water brake dynamometer. In the past this 
inertia has been assumed to be a constant value when, in fact, it 
actually changes with the addition and subtraction of water. 

The last recommendations for this project encompass the entire 
project. First, a digital data acquisition system should be used to get 
the steady state and dynamic data. This will make the fine tuning 
process with the data regression program more accurate. Finally, one 
computer system should be used for the whole process, from the 
regression and grid search to the computer simulation (vice the two 
system used in this project). In this way, results could be achieved 
and correlated in a more efficient manner. 

The model was robust in that, with the exception of the constant 
Ns in the vertical runs, the model was applicable to all areas of the 
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operating range. This was supported by the excellent response in Ng 
throughout the whole operating range. This shows that the 

coefficients for the All, Al2, and A13 were accurate and the 
procedure to find them was sound. 

Once the proper values for A21, A22, and A23 are found and this 

modeling technique has been successfully demonstrated on the Boeing 
engine this method will have real time application to the future 
controller design for the General Electric LM- 2500 and the follow on 
marine gas turbine controllers designed for the future. 
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APPENDIX A 



COMPONENT MODELING EQUATIONS : 



Compressor 




dt2 = Al dng + A2 dp2 
dma = Bl dng + B2 dp2 
dqc = Cl dng + C2 dp2 



1 

2 

3 



Gas Generator Shaft 




qhpt JGB dng = dqhpt-dqc 4 

ng 



Dynamometer Shaft 




JDB dns = dqfpt-dqd 



5 
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Free Power Turbine 




dp4 = Dl dng + D2 dns 6 



dqfpt = Fl dns+F2 dt4+F3 dmaf ... 7 



Hi gh Pressure Turbine 




dt4 = Gl dma + G2 dt2 + G3 dp4 + G4 dng + G5 de 10 

dp2 = Hl dng + H2 dns 11 

dqhpt = U dma + 12 dt2 + 13 dp4 + 14 dng + 15 de 12 





Note: 1. Lower case variables represent normalized variables. 
2. The d _ represents the perturbation of the variable. 

• Decoupling equations. 



42 



APPENDIX B 



EQUATION SET REDUCTION METHODOLOGY: 

The equations in appendix A are reduced into the state space 
components of the A and B matrices. 



Substitute eqn 1 1 : 

dp2 = Hi dng + H2 dns 1 1 

into eqns 1, 2, and 3 for dp2. 

dt2 = (Al + A2Hl) dng + (A2H2) dns 1 

dma = (Bl + B2H1) dng + (B2H2) dns 2 

dqc = (Cl + C2Hl) dng + (C2H2) dns 3 

In eqn 12: 

dqhpt = II dma + 12 dt2 + 13 dp4 + 14 dng + 15 de 12 

substitute eqns 1, 2, and 6 for dt2, dma, and dp4. 
dqhpt = 11 (Bl + B2H1) dng + U (B2H2) dns + 

12 (Al + A2H1) dng + 12 (A2H2) dns + 

13 (Dl) dng + 13 (D2) dns + 

14 dng + 15 de 12 
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Collect terms: 



dqhpt = (U (Bl + B2H1) + 12 (Al + A2H1) + 13 (Dl) + 14 ) dng + 
(11 (B2H2) + 12 (A2H2) + 13 (D2) ) dns + 15 de 

Let Y1 = II (Bl + B2H1) + 12 (Al + A2Hl) + 13 (Dl) + 14 
Y2 = II (B2H2) + 12 (A2H2) + 13 (D2) 

Then: dqhpt = Y1 dng + Y2 dns + 15 de 12 

Take eqn 4 and substitute in new eqn 12 for dqhpt and eqn 3 for 
dqc. 

JGB dng = Y1 dng + Y2 dns + 15 de - 

(Cl + C2Hl) dng - (C2H2) dns 
Collect terms: 

dng = ( (Y1 - (Cl + C2H1)) dng + (Y2 - (C2H2)) dns + 15 de )/JGB 

Let: Z1 = (Y1 - (Cl + C2H1) 

Z2 = Y2 - (C2H2) 

Z3 = 15 

Then: 

dng = (Z1 dng + Z2 dns + Z3 de)/JGB R1 
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In eqn 10: 

dt4 = Gl dma + G2 dt2 + G3 dp4 + G4 dng + G5 de 10 

substitute eqns 1, 2, and 6 for dt2, dma, and dp4. 
dt4 = Gl (B1 + B2H1) dng + Gl (B2H2) dns + 

G2 (Al + A2H1) dng + G2 (A2H2) dns + 

G3 Pi) dng + G3 (D2) dns + 

G4 dng + G5 de 10 

Collect terms: 

dt4 = (Gl (Bl + B2H1) + G2 (Al + A2Hl) + G3 (Dl) + G4 ) dng + 
(Gl (B2H2) + G2 (A2H2) + G3 (D2) ) dns + G5 de 

Let Y3 = Gl (Bl + B2Hl) + G2 (Al + A2Hl) + G3 (Dl) + G4 
Y4 = Gl (B2H2) + G2 (A2H2) + G3 (D2) 

Then: dt4 = Y3 dng + Y4 dns + G5 de 10 

In the continuity eqn (9) substitute eqn 2 for dma. 

dmaf = Kl dmf + K2(B1 + B2H1) dng + K2(B2H2) dns 9 

Substitute new eqn 10 for dt4 and new eqn 9 for dmaf into: 
eqn 7: dqfpt = Fl dns + F2 dt4 + F3 dmaf 

This results in: 
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dqfpt = Fl dns + F2Y3 dng + F2Y4 dns + F2G5 de + 

F 3 K 1 dmf + F 3 K 2 (B 1 + B 2 H 1 ) dng + F3K2(B2H2) dns 
Collect terms: 

dqfpt = (F2Y3 + F 3 K 2 (Bl + B 2 H 1 )) dng + 

(Fl + F3K2(B2H2) + F2Y4) dns + F2G5 de + F 3 K 1 dmf 

Let: Z4 = F2Y3 + F 3 K 2 (Bl + B 2 H 1 ) 

Z5 = Fl + F3K2(B2H2) + F2Y4 
Z 6 = F2G5 
Z7 = F 3 K 1 
Then: 

dqfpt = Z4 dng + Z5 dns + Z 6 de + Z7 dmf 7 

Substitute new eqn 7 for qfpt in eqn 5. 



JDB dns = dqfpt - dqd 5 

dns = (Z4 dng + Z5 dns + Z6 de + Z7 dmf - dqd) /JDB R2 



Finally: Eqn 8 is R3. Rl, R2, and R3 make up the state space 
equation set. 
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The state space equation set is the following: 



dng = Zi dng + Z2 dns + Z2 de R1 

JGB JGB JGB 

dns = Z4 dng + Z£ dns + ZQ. de + £7 dmf - 1 dqd R2 

JDB JDB JDB JDB JDB 

de = - 1 de + 1 dmf R3 

T T 



The resulting A and B matrices are: 


All =Z1 


A12 = Z2 


A13 = Z3 


JGB 


JGB 


JGB 


A21 = Z4 


A22 = Z5 


A23 = Z6 


JDB 


JDB 


JDB 


A31 = O 


A32 = O 


A33 = J_ 



T 



and 

Bll = O B12 = O 



B21=Z7 B22 = -1 

JDB JDB 



B31 = JL B32 = O 
T 



4 7 



APPENDIX £ 



PARAMETER SELECTION METHODOLOGY: 

The parameter selections were accomplished through three means: 

1 . Plotting actual slopes from data. This was shown 
in chapter III section F. 

2. Data reduction using mainframe program 
Mini-Tab. 

3. Grid searching using coefficients from 1. and 2. 
and known relationships which must be satisfied. 
The spread sheet is from the Micro-soft program 

Excell and was used on an Apple Macintosh 
computer. 

The data reduction program Minitab was used for all of the 
regressions performed. The input for the regression program was the 
steady state turbine data shown in table 1 A. This routine was used to 
determine the order of the independent variables in the equations for 
Ma, T2, and Qc, as well as to regress the expressions for P2, P4, T4, 
Qhpt and Maf. To determine the order of a independent variable a 
three step process was followed. First, the dependent variable was 
regressed using the data for the independent variables. In doing this a 
residual was calculated for each set of the independent variables. This 
was then subtracted from the regressed value of the dependent 
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variable. The resulting difference was plotted versus each of the 
independent variables. In the plot for the regression of T2, shown in 
figure C-la, there is a definite parabolic shape which is an indication 
that the independent variable in question (Ng) is of a higher order. To 

verify this, T2 was regressed again, this time using Ng A ^, Ng and P2. 
The plotting was repeated with the resultant plot shown in figure C- 
lb. This time the parabolic shape was not in evidence and the 
predictor percentage, which indicates the accuracy of the regression 
based on the data, increased from 99.6 to 99.9 percent indicating an 
excellent regression for the t2 equation. The resulting partials for all 
of the component expressions are shown in figure C-2 lines 3 through 
21 . 

The spreadsheet grid, figure C-2, was used in conjunction with the 
simulation program to search out the best values for the rest of the 
coefficients in the equations in appendix A. These were input for Gl, 
F2, and F3 in lines 3-21 and then used in lines 25-32 to calculate the 
values of the coefficients of the A and B matrices. Lines 34-38 show 
the proper signs for these matrix coefficients as well as the dns and 
dng relationships which must be satisfied during the search process. 
Once these were satisfied, plots were generated and checked for the 
accuracy of the simulation. 
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Figure C-l. Regression Plots 
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Figure C-2. Sample Grid Search Spread Sheet. 



APPENDIX D 



COMPUTER SIMULATION CODES: 



The computer codes were written using the Dynamic Simulation 
Language (DSL). 

Simulation Run 1: DSL code for dynamic Ns run at constant Ng 

(20000 rpm). 

Simulation Run 3: DSL code for dynamic Ng run at constant Ns 

(500 rpm). 

Simulation Run 5: DSL code for dynamic Ng run at constant Ns 

(1500 rpm). 

Simulation Run 7 : DSL code for dynamic Ng run at Ns 

( 2252-2500 rpm). 

Simulation Run 9: DSL code for dynamic Ns run at constant Ng 

(32000 rpm). 
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********************************************************************** 



* 



★ 



★ 



★ 



* BOEING MODEL 502-6A GAS TURBINE 

* 



* 

* 



* 

★ 



DYNAMIC COMPUTER SIMULATION * 

* 



* MODIFIED BY 

* S.D. METZ 

* FROM ROUTINES BY 

* V. J.HERDA AND V. A. STAMMETTI 



* 

★ 

* 

* 

* 



* THIS PROGRAM SIMULATES THE DYNAMIC RESPONSE OF THE NPS * 

* BOEING GAS TURBINE TEST FACILITY USING A MULTI LPLE * 

* LINEARIZATION TECHNIQUE. * 

* * 



*********************************************************************** 



C 

c 

PARAM JG=0. 009525, JD=0.3000, PI=3. 14159, TC=1.0 

* 



* THE FOLLOWING VALUES LISTED ON THE FUNCTION CARD ARE FOR FUEL FLOW, 

* GAS GENERATOR SPEED, TORQUE AND DYNO SPEED AS A FUNCTION OF TIME. 

* THESE VALUES WERE OBTAINED FROM STRIP CHART RECORDS AND ARE ENTERED 

* IN THE FORM ( E. G. FUEL FLOW) . . . TIME( SEC ) , FUEL FLOW 



C 

C THIS SET IS FOR EXPERIMENTAL RUN # 1. 

C 

C 

AFGEN NGDATA= 0.0,20230.0, 5.0,20230.0, 10.0,20230.0, 15.0,20230.0, ... 

20. 0. 20230. 0, 25. 0,20230. 0, 50. 0,20230. 0 

AFGEN NSDATA= 0. 0, 505. 0, 1.0,505.0, 2.0,536.7, 3.0,568.33, ... 

4.0. 600.0, 5.0,631.7, 6.0,695.0, 7.0,790.0, 8.0,885.0, ... 

9.0. 1043.3, 10.0,1170.0, 11.0,1328.3, 12.0,1486.7, ... 

13.0. 1708.3, 14.0,1835.0, 15.0,1930.0, 16.0,1993.3, ... 

17.0. 2056.7, 18.0,2088.3, 19.0,2120.0, 20.0,2151.7, ... 

21.0. 2167.5, 22.0,2177.0, 23.0,2183.3, 24.0,2215.0, ... 

25. 0. 2215. 0, 50. 0,2215. 0 

AFGEN MFDATA= 0. 0,78. 39, 5.0,78.39, 10.0,78.39, 12.0,79.89, ... 

15.0. 80.26, 17.0,80.26, 19.0,81.38, 20.0,82.50, ... 

21.0. 82.13, 25.0,82.13, 50.0,82.13 



AFGEN QDDATA= 0. 0, 125. 0, 1.0,125.0, 2.0,118.3, 3.0,115.0, ... 

4.0. 111.7, 5.0,105.0, 6.0,98.33, 7.0,85.0, ... 

8.0. 71.67, 9.0,58.33, 10.0,45.0, 11.0,31.67, ... 

12.0. 11.67, 12.5,1.670, 13.0,5.0, 15.0,5.0, ... 

25.0. 5.0, 50.0,5.0 



* 

INITIAL 



* 

* 



* 



* 

* 

* 



ESTABLISH INITIAL CONDITIONS. 



SET 



NGI =20230. 0 
NSI =505. 0 



MFI 


= 


78. 39 


QDI 




125. 0 


INITIAL 


STATE 


DNG 


= 


0. 0 


DNS 


= 


0. 0 


DE 


= 


0. 0 



PERTURBATION TO ZERO 
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NGB = 2 6000.0 

NSB = 1500. 0 

MFB = 122. 5 

QDB = 220.0 

QCB = 80. 0 

JGB = ( NGB/QCB ) * JG 

JDB = ( NSB/QDB ) * JD 



DYNAMIC 

★ 



DATA CURVE FORMULATION 

NGD = AFGEN( NGDATA, TIME ) 

NSD = AFGEN( NSDATA, TIME ) 

MFD = AFGEN( MFDATA , TIME ) 

QDD = AFGEN( QDDATA, TIME) 

MFDLY = TRANSP( 100 , MFI , 0. 70, MFD ) 

STATE SPACE LINEAR MODEL FORMULATION 
All = < 20. 282-32. 560*( NGL/NGB) )/( JGB) 
A12 = ( -0. 4480 )/( JGB) 

A13 = ( 13. 624 )/( JGB) 

A21 = ( 2. 1131 + . 91786*( NGL/NGB) )/( JDB) 
A22 = (-. 0998-0. 445*( NGL/NGB) )/( JDB) 
A23 = { -0. 1885 ) /( JDB) 

A3 1 = 0. 0 

A32 = 0. 0 

A33 = -1. O/TC 

Bll = 0. 0 

B12 = 0. 0 

B21 =. 035250/( JDB) 

B22 = -1. 0/( JDB) 

B31 = 1. O/TC 
B32 = 0. 0 



DERIVATIVE 

NOSORT 

* 

* COMPUTE INPUT TO THE NONLINEAR MODEL, MF( T ) , WW( T ) . 

* RUN #9 

DMF =( MFDLY - MFI) /MFB 
DQD =( QDD - QDI)/QDB 

DNGDOT = All *DNG + A12*DNS + A13*DE + B11*DMF 
DNSDOT = A2 1 *DNG + A22*DNS + A23*DE + B21*DMF ■« 
DEDOT = A33 *DE + B3 1 *DMF 



DYNAMIC EQUATIONS FOR LINEAR MODEL. 
DNG= I NTGRL( 0.0, DNGDOT ) 
DNS=INTGRL( 0. 0, DNSDOT) 

DE =INTGRL( 0. 0, DEDOT) 

NGL = NGI + DNG^NGB 
NSL = NS I + DNS*NSB 



SORT 



B22*DQD 
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* THE STATEMENTS IN THE PREVIOUS (DERIVATIVE) SECTION YIELD VALUES 

* OF ' NG ' , AND ' NS 1 AS CALCULATED BY THE NONLINEAR AND STATE SPACE 

* MODELS. THE STATEMENTS BELOW COMPUTE THE VALUES OF 'NG', AND 'NS' 

* AS RECORDED FROM GAS TURBINE TEST DATA. 

* 

TERMINAL 
METHOD RK5 

CONTROL FINTIM=50. 0,DELT=0. 001 
PRINT 0. 5 , NGD, NGL, NSD, NSL, MFD, QDD 
END 
STOP 
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*+****★**+★*★★★*★★*★*★★★**★***++*+***+++*++**+*★**★★*******+*++****«*« 



★ 



★ 



★ 

* 

★ 

★ 



BOEING MODEL 502-6A GAS TURBINE * 

* 

DYNAMIC COMPUTER SIMULATION * 

★ 



* MODIFIED BY * 

* S. D. METZ * 

* FROM ROUTINES BY * 

* V. J.HERDA AND V. A. STAMMETTI * 

* * 

* THIS PROGRAM SIMULATES THE DYNAMIC RESPONSE OF THE NPS * 

* BOEING GAS TURBINE TEST FACILITY USING A MULTILPLE * 

* LINEARIZATION TECHNIQUE. * 

* * 



it********************************************************************** 



C 

c 

PARAM JG=0. 009525, JD=0. 3000, PI=3. 14159, TC=1. 00 

* 

* THE FOLLOWING VALUES LISTED ON THE FUNCTION CARD ARE FOR FUEL FLOW, 

* GAS GENERATOR SPEED, TORQUE AND DYNO SPEED AS A FUNCTION OF TIME. 

* THESE VALUES WERE OBTAINED FROM STRIP CHART RECORDS AND ARE ENTERED 

* IN THE FORM ( E. G. FUEL FLOW) . . . TIME( SEC ) , FUEL FLOW 



C 

C THIS SET IS FOR EXPERIMENTAL RUN # 3. 

C 

C 

AFGEN NGDATA= 0. 0,20220. 0, 3.55,20220.0, 3.7,20315.7, 3.8,20411.4, ... 

4.0,20698.4, 5.0,23186.1, 6.0,25482.4, 7.0,27204.6, ... 

8.0,28926.9, 9.0,30649.1, 10.0,32180.0, 10.5,32371.4, ... 
11. 0,32180. 0, 11. 5,3 2084. 3, 12. 0,32180. 0, 50. 0,32180. 0 

AFGEN NSDATA= 0. 0,506. 0, 2.0,505.7, 3.6,505.5, 4.0,506.5, ... 

5.0. 508.3, 6.0,509.7, 7.0,511.1, 8.0,512.5, 9.0,513.5, ... 

10.0. 514.4, 10.4,514.7, 11.0,514.4, 12.0,515.1, ... 

13.0. 515.8, 14.0,516.7, 16.3,517.2, 17.0,516.7, ... 

18.0. 515.8, 19.0,514.9, 20.0,513.0, 50.0,513.0 
AFGEN MFDATA= 0. 0, 77. 5, 3.0,77.5, 3.3,79.6, 3.6,85.8, ... 

4.0. 91.99, 5.0,106.5, 6.0,118.9, 7.0,130.5, 8.0,141.3, ... 

9.0. 156.2, 10.0,172.7, 11.0,156.2, 12.0,160.3, 50.0,160.3 



AFGEN QDDATA= 0. 0, 125. 0, 3.6,125.0, 3.8,128.4, 4.0,131.8, ... 

5.0. 172.4, 6.0,226.7, 7.0,274.1, 8.0,321.6, 9.0,375.8, ... 

10.0. 430.0, 10.5,443.6, 11.0,430.0, 12.0,416.4, ... 

13.0. 416.4, 14.0,419.8, 15.0,416.4, 18.0,423.0, ... 

20. 0. 430. 0, 50. 0, 430 



INITIAL 

* ESTABLISH INITIAL CONDITIONS. 

* 

NGI=20220. 0 
NSI=506. 0 
MFI = 77. 5 
QDI = 125. 0 



* SET INITIAL STATE PERTURBATION TO ZERO 

★ 

DNG = 0.0 
DNS = 0.0 
DE = 0. 0 
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NGB = 26000. 0 

NSB = 1500.0 

MFB = 122. 5 

QDB = 220. 0 

QCB = 80. 0 

JGB = ( NGB/QCB) * JG 

JDB = ( NSB/QDB ) * JD 

DYNAMIC 

★ 

* DATA CURVE FORMULATION 

* 

NGD = AFGEN(NGDATA,TIME) 

NSD = AFGEN( NSDATA, TIME) 

MFD = AFGEN( MFDATA, TIME ) 

QDD = AFGEN( QDDATA, TIME ) 

MFDLY = TRANSP( 100, MF 1 , 0. 70, MFD ) 



* 

* 



* 

★ 

★ 



* 



* 



STATE SPACE LINEAR MODEL FORMULATION 
All = ( 20. 282-32. 560*( NGL/NGB) )/( JGB) 
A12 = ( -0. 448 )/( JGB) 

A13 = ( 13. 624 )/( JGB) 

A21 = (2. 1131+. 91786*( NGL/NGB) )/( JDB) 
A22 = (-. 0998-0. 445*( NGL/NGB) )/( JDB) 
A23 =( -. 18851 )/( JDB) 

A3 1 = 0. 0 

A32 = 0.0 

A3 3 = -1.0/TC 

Bll = 0. 0 

B12 = 0. 0 

B21 =. 035250/( JDB) 

B22 = -1. 0/( JDB) 

B31 = 1.0/TC 
B32 = 0. 0 



DERIVATIVE 

NOSORT 

* 

* COMPUTE INPUT TO THE NONLINEAR MODEL, MF( T ) , WW(T). 

* RUN #3 

DMF =( MFDLY - MFI ) /MFB 
DQD =( QDD - QDI)/QDB 

DNGDOT = All *DNG + A12*DNS + A13*DE + Bll+DMF 
DNSDOT = A2 1 *DNG + A22*DNS + A23*DE + B2 1 *DMF 
DEDOT = A33 *DE + B3 1 *DMF 



* DYNAMIC EQUATIONS FOR LINEAR MODEL. 
DNG=INTGRL( 0. 0, DNGDOT) 

DNS = INTGRL( 0. 0, DNSDOT) 

DE =INTGRL( 0.0, DEDOT ) 

NGL = NGI + DNG*NGB 
NSL = NS I + DNS*NSB 

SORT 



B22 *DQD 



57 



* THE STATEMENTS IN THE PREVIOUS (DERIVATIVE) SECTION YIELD VALUES 

* OF 'NG', AND 'NS' AS CALCULATED BY THE NONLINEAR AND STATE SPACE 

* MODELS. THE STATEMENTS BELOW COMPUTE THE VALUES OF 'NG' , AND 'NS' 

* AS RECORDED FROM GAS TURBINE TEST DATA. 

* 

TERMINAL 
METHOD RK5 

CONTROL FINTIM=50. 0,DELT=0. 001 

PRINT 0. 50,NGD,NGL,NSD,NSL,MFD,QDD,DNSD0T,A21,DNG, . . . 

A22 , DNS , A23 , DE , B21 , DMF, B22 , DQD 

END 

STOP 
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+**+*+**++★*++****★*++**+***★+*****★*++*+++*+++********+*+*+***+**★★★* 
* * 



* BOEING MODEL 502-6A GAS TURBINE 

* 

* DYNAMIC COMPUTER SIMULATION 

* 

* MODIFIED BY 

* S.D. METZ 

* FROM ROUTINES BY 

* V. J.HERDA AND V. A. STAMMETTI 

* 



* 

* 



★ 

* 



* THIS PROGRAM SIMULATES THE DYNAMIC RESPONSE OF THE NPS * 

* BOEING GAS TURBINE TEST FACILITY USING A MULTILPLE * 

* LINEARIZATION TECHNIQUE. * 

* * 



*********************************************************************** 



C 

C 

PARAM JG=0. 009525, JD=0. 3000, PI=3. 14159, TC=1. 00 
★ 

* THE FOLLOWING VALUES LISTED ON THE FUNCTION CARD ARE FOR FUEL FLOW, 

* GAS GENERATOR SPEED, TORQUE AND DYNO SPEED AS A FUNCTION OF TIME. 

* THESE VALUES WERE OBTAINED FROM STRIP CHART RECORDS AND ARE ENTERED 

* IN THE FORM ( E. G. FUEL FLOW) . . . TIME( SEC ) , FUEL FLOW 



C 

C THIS SET IS FOR EXPERIMENTAL RUN # 5. 

C 

C 

AFGEN NGDATA= 0. 0,20220. 0, 3.1,20220.0, 3.6,20412.1, 4.0,21372.6, ... 

5.0. 23677.7, 6.0,25886.9, 7.0,27711.8, 8.0,29440.6, ... 

9.0. 30977.4, 10.0,32206.8, 10.2,32322.1, 10.8,32130.0, ... 

11.5,32130.0, 12.1,32014.7, 13.4,32206.8, 14.5,32014.7, ... 
15. 6,32206. 8, 16. 8, 32014. 7 , 17. 0,32206. 8, 50. 0,32130. 0 

★ 

AFGEN NSDATA= 0. 0, 1491. 6, 1.0,1490.4, 3.0,1489.6, 4.0,1492.5, ... 

5.0. 1501.0, 6.0,1508.3, 7.0, 1511.7, 7.8,1513.3, ... 
9.6,1513.3, 10.1,1513.8, 12.0,1508.3, 13.0,1510.0, ... 

50. 0. 1510. 0 

* 

AFGEN MFDATA= 0. 0,83. 75, 3.2,83.75, 5.0,109.6, 6.0,122.2, ... 

7.0. 135.5, 8.0,150.2, 9.0,162.8, 9.8,172.4, ... 

10.0. 168.7, 11.0,165.0, 50.0,165.0 

★ 

AFGEN QDDATA= 0. 0,71. 0, 4.4,71.0, 5.0,78.7, 6.0,94.1, ... 

7.0. 107.5, 8.0,120.9, 9.0,134.4, 10.2,149.8, ... 

10.9,147.9, 11.6,147.9, 12.4,143.2, 13.0,145.9, ... 

13.5.144.0, 14.0,145.9, 15.0,144.0, 50.0,144.0 



INITIAL 

* ESTABLISH INITIAL CONDITIONS 

* 



★ 

★ 

* 

* 

★ 



SET 



NGI=20220. 0 
NSI=1491. 6 



MFI 


= 


83. 75 


QDI 


— 


71. 0 


INITIAL 


STATE 


DNG 


= 


0. 0 


DNS 


= 


0. 0 


DE 


= 


0. 0 



PERTURBATION TO ZERO 
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DYNAMIC 

★ 



NGB = 26000.0 

NSB = 1500. 0 

MFB = 122. 5 

QDB = 220.0 

QCB = 60. 0 

JGB = ( NGB/QCB) * JG 

JDB = ( NSB/QDB ) * JD 



DATA CURVE FORMULATION 



NGD = AFGEN( NGDATA, TIME ) 

NSD = AFGEN( NSDATA, TIME ) 

MFD = AFGEN( MFDATA , TIME ) 

QDD = AFGEN( QDDATA, TIME ) 

MFDLY = TRANSP( 100,MFI,0. 70, MFD) 

STATE SPACE LINEAR MODEL FORMULATION 

All ■ (20. 282-32. 560*( NGL/NGB) )/( JGB) 



A12 


=2 


( -0. 4480 ) /( JGB) 


A13 


= 


( 13. 624)/( JGB) 


A21 


= 


(2. 1 13 1 + . 91786*( NGL/NGB) )/( JDB) 


A22 


= 


( 0998-0. 445*( NGL/NGB) )/( JDB) 


A23 


= ( 


-0. 18851 )/( JDB) 


A3 1 


= 


0. 0 


A32 


= 


0. 0 


A33 


= 


-1. 0/TC 


Bll 


= 


0. 0 


B12 


= 


0. 0 


B21 


= . 


03 5250/( JDB) 


B22 


= 


-1. 0/( JDB) 


B31 


= 


1. 0/TC 


B32 


= 


0. 0 



DERIVATIVE 

NOSORT 

★ 

* COMPUTE INPUT TO THE NONLINEAR MODEL, MF( T ) , WW(T). 

* RUN #9 

DMF =( MFDLY - MFI ) /MFB 
DQD =( QDD - QDI)/QDB 

DNGDOT = A1 1 *DNG + A12*DNS + A13*DE + B11*DMF 

DNSDOT = A2 1 *DNG + A22*DNS + A23*DE + B21*DMF + B22*DQD 

DEDOT = A33 *DE + B3 1 *DMF 



DYNAMIC EQUATIONS FOR LINEAR MODEL. 
DNG=INTGRL( 0. 0, DNGDOT) 

DNS= INTGRL( 0. 0, DNSDOT) 

DE =INTGRL( 0.0, DEDOT ) 

NGL = NGI + DNG*MGB 
NSL = NS I + DNS*NSB 
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* THE STATEMENTS IN THE PREVIOUS (DERIVATIVE) SECTION YIELD VALUES 

* OF 'NG', AND ’NS' AS CALCULATED BY THE NONLINEAR AND STATE SPACE 

* MODELS. THE STATEMENTS BELOW COMPUTE THE VALUES OF ' NG ' , AND 'NS' 

* AS RECORDED FROM CAS TURBINE TEST DATA. 

* 

TERMINAL 
METHOD RK5 

CONTROL FINTIM=50. 0,DELT=0. 001 

PRINT 0.5, NGD , NGL , NSD , NSL , MFD , QDD , DNSDOT , A2 1 , DNG , A2 2 , DNS , ... 

A2 3 , DE , B2 1 , DMF , B22 , DQD 

END 

STOP 
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********************************************************************** 
* * 



* * 

* BOEING MODEL 502-6A GAS TURBINE * 

★ * 

* DYNAMIC COMPUTER SIMULATION * 

★ ★ 

* MODIFIED BY * 

* S.D. METZ * 

* FROM ROUTINES BY * 

* V. J.HERDA AND V. A. STAMMETTI * 

* * 



* THIS PROGRAM SIMULATES THE DYNAMIC RESPONSE OF THE NPS * 

* BOEING GAS TURBINE TEST FACILITY USING A MULTILPLE * 

* LINEARIZATION TECHNIQUE. * 

* ★ 



★*★★★*★★*★★★★*★★*★★★★★★★★★*★★★*★★★★*★*★*★ ★★★★*★★★***★★★★****★***★★**★*★ 



C 

C 

PARAM JG=0. 009525, JD=0.3000, PI=3. 14159, TC=1. 00 
★ 



* THE FOLLOWING VALUES LISTED ON THE FUNCTION CARD ARE FOR FUEL FLOW, 

* GAS GENERATOR SPEED, TORQUE AND DYNO SPEED AS A FUNCTION OF TIME. 

* THESE VALUES WERE OBTAINED FROM STRIP CHART RECORDS AND ARE ENTERED 

* IN THE FORM ( E. G. FUEL FLOW) . . . TIME( SEC ) , FUEL FLOW 



C 

C THIS SET IS FOR EXPERIMENTAL RUN # 7. 

C 

C 

AFGEN NGDATA= 0. 0,20220. 0, 4.0,20220.0, 4.5,20759.7, 5.0,21954.7, ... 

6.0,24460.3, 7.0,26541.9, 8.0,28122.4, 9.0,29780.0, ... 

10.0. 31399.0, 10.4,32092.9, 10.6,32362.7, 11.0,32285.6, ... 

12. 0. 32170. 0, 15. 0,32170. 0, 20. 0,32170 

AFGEN NSDATA= 0. 0,2252. 0, 1.0,2245.4, 2.0,2238.8, 3.0,2232.2, ... 

4.0. 2225.6, 5.0,2278.4, 6.0,2436.8, 7.0,2568.8, ... 
7.7,2621.6, 8.0,2611.0, 9.0,2579.4, 10.0,2568.8, ... 

11.0. 2542.4, 13.0,2516.0, 15.0,2516.0, 20.0,2516 
AFGEN MFDATA= 0. 0, 78. 8, 1.0,78.8, 3.0,78.8, 4.0,87.1, ... 

4.5,95.3, 5.0,103.6, 6.0,111.9, 7.0,128.5, 8.0,136.8, ... 

9.0. 161.7, 10.0,174.1, 10.4,178.3, 11.0,170.0, ... 

15.0. 170.0, 20.0,170.0 

AFGEN QDDATA= 0. 0,2 5. 0, 2.0,25.0, 4.0,25.0, 5.0,25.0, ... 

6.0. 49.6, 7.0,92.6, 8.0,129.4, 9.0,178.6, ... 

10.0. 233.9, 11.0,264.6, 12.0,252.3, 13.0,246.1, ... 

15.0. 240.0, 20.0,240.0 



INITIAL 

★ 

* ESTABLISH INITIAL CONDITIONS. 

NGI=20220. 0 
NSI=2252. 0 
MFI = 78. 75 
QDI = 25. 0 

★ 

★ 

★ 

* SET INITIAL STATE FERTURBATION TO ZERO 

★ 

DNG = 0.0 
DNS = 0. 0 

DE = 0. 0 
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DYNAMIC 

★ 



NGB = 26000.0 

NSB = 1500.0 

MFB = 122. 5 

QDB = 220. 0 

QCB = 80. 0 

JGB = ( NGB/QCB) * JG 

JDB = ( NSB/QDB ) * JD 



DATA CURVE FORMULATION 



NGD = AFGEN( NGDATA , T I ME ) 

NSD = AFGEN( NSDATA, TIME) 

MFD = AFGEN( MFDATA, TIME ) 

QDD = AFGEN( QDDATA, TIME) 

MFDLY = TRANSP( 100,MFI , 0. 70, MFD) 

STATE SPACE LINEAR MODEL FORMULATION 
All = (20. 282-32. 560* ( NGL/NGB) )/( JGB) 
A12 = (-0. 448)/( JGB) 

A13 = ( 13. 624)/( JGB) 

A21 = ( 2. 1131+. 91786*( NGL/NGB) )/( JDB) 
A22 = (-. 0998-0. 445*( NGL/NGB) )/( JDB) 
A23 = ( -. 18851 )/( JDB) 

A3 1 =0.0 

A32 = 0. 0 

A33 = -1. O/TC 

Bll = 0. 0 

B12 = 0. 0 

B21 =. 035250/( JDB) 

B22 = -1. 0/( JDB) 

B31 = 1. O/TC 
B32 =0.0 



* 

* 

DERIVATIVE 

NOSORT 

* 

* 

* 

* 



COMPUTE INPUT TO THE NONLINEAR MODEL, MF(T), WW(T). 

RUN #7 

DMF =( MFDLY - MFI )/MFB 
DQD =( QDD - QDI ) /QDB 

DNGDOT = A11*DNG + A12*DNS + A13*DE + B11*DMF 

DNSDOT = A21*DNG + A22*DNS + A23*DE + B21*DMF + B22*DQD 

DEDOT = A33 *DE + B31 *DMF 



DYNAMIC EQUATIONS FOR LINEAR MODEL. 
DNG=INTGRL( 0. 0, DNGDOT) 
DNS=INTGRL( 0. 0, DNSDOT) 

DE =INTGRL( 0. 0, DEDOT) 

NGL = MG I + DNG*NGB 
NSL = NS I + DNS* NSB 



SORT 

* 
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* OF 'MG', AND 'NS' AS CALCULATED BY THE NONLINEAR AND STATE SPACE 

* MODELS. THE STATEMENTS BELOW COMPUTE THE VALUES OF 'NG', AND ’NS’ 

* AS RECORDED FROM GAS TURBINE TEST DATA. 

* 

TERMINAL 
METHOD RK5 

CONTROL FINTIM=20. 0,DELT=0. 001 
PRINT 0. 25,NGD,NGL,NSD,NSL,MFD,QDD 
END 
STOP 
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ft********************************************************************* 



* * 

* * 

* BOEING MODEL 502-6A GAS TURBINE * 

* * 

* DYNAMIC COMPUTER SIMULATION * 

★ * 

* MODIFIED BY * 

* S.D. METZ * 

* FROM ROUTINES BY * 

* V. J.HERDA AND V. A. STAMMETTI * 

* * 



* THIS PROGRAM SIMULATES THE DYNAMIC RESPONSE OF THE NPS * 

* BOEING GAS TURBINE TEST FACILITY USING A MULTILPLE * 

* LINEARIZATION TECHNIQUE. * 

* * 
*********************************************************************** 
C 

C 

PARAM JG=0. 009525, JD=0. 3000, PI=3. 14159, TC=1. 00 
* 

* THE FOLLOWING VALUES LISTED ON THE FUNCTION CARD ARE FOR FUEL FLOW, 

* GAS GENERATOR SPEED, TORQUE AND DYNO SPEED AS A FUNCTION OF TIME. 

* THESE VALUES WERE OBTAINED FROM STRIP CHART RECORDS AND ARE ENTERED 

* IN THE FORM ( E. G. FUEL FLOW) . . . TIME( SEC ) , FUEL FLOW 



C 

C THIS SET IS FOR EXPERIMENTAL RUN # 9. 

C 

c 

AFGEN NGDATA= 0. 0,32150. 0, 5.0,32150.0, 10.0,32150.0, 15.0,32150.0, ... 

20.0. 32150.0, 25.0,32150.0, 50.0,32150.0 

AFGEN NSDATA= 0. 0, 510. 0, 2.0,510.0, 3.0,576.1, 4.0,642.14, ... 

5.0. 708.2, 6.0,774.3, 7.0,840.4, 8.0,922.9, 9.0,1005.5, ... 

10.0. 1137.7, 11.0,1269.8, 12.0,1435.0, 13.0,1583.6, ... 

14.0. 1798.4, 15.0,1963.5, 16.0,2128.7, 17.0,2227.8, ... 

18.0. 2326.9, 19.0,2369.0, 20.0,2409.5, 21.0,2459.1, ... 

22.0. 2492.1, 23.0,2508.6, 24.0,2525.0, 25.0,2525.0, ... 

50. 0. 2525. 0 

AFGEN MFDATA= 0. 0, 155. 8, 3.0,155.8, 5.0,155.8, 7.0,160.0, ... 

10.0. 160.0, 13.0,164.1, 15.0,164.1, 18.0,164.1, ... 

20.0. 168.3, 23.0,168.3, 24.0,168.3, 25.0,168.3, 50.0,168.3 



AFGEN QDDATA= 0. 0, 394. 0, 2.0,394.0, 3.0,382.5, 4.0,371.0, ... 

5.0. 359.5, 7.0,336.5, 8.0,325.0, 9.0,313.50, ... 

10.0. 302.0, 11.0,284.8, 12.0,261.8, 13.0,250.3, ... 

14.0. 233.0, 15.0,221.5, 16.0,221.5, 19.0,221.5, ... 

20.0. 221.5, 21.0,210.0, 23.0,210.0, 25.0,210.0, 50.0,210.0 



INITIAL 



* ESTABLISH INITIAL CONDITIONS. 

* 

NGI=32 150. 0 
NS 1 = 510. 0 
MFI = 155. 8 
QDI = 394. 0 



* SET INITIAL STATE PERTURBATION TO ZERO 

* 

DNG = 0.0 
DNS = 0.0 
DE =0.0 
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NGB = 26000.0 

NSB = 1500. 0 

MFB = 122. 5 

QDB = 220. 0 

QCB = 80. 0 

JGB = ( NGB/QCB ) * JG 

JDB = ( NSB/QDB ) * JD 

DYNAMIC 

★ 

* DATA CURVE FORMULATION 

* 

NGD = AFGEN( NGDATA, TIME ) 

NSD = AFGEN( NSDATA, TIME ) 

MFD = AFGEN( MFDATA, TIME ) 

QDD = AFGEN( QDDATA, TIME ) 

MFDLY = TRANSP( 100,MFI,0. 70, MFD) 

* 

* STATE SPACE LINEAR MODEL FORMULATION 

* 

All = ( 20. 282-32. 560* ( NGL/NGB ) )/( JGB) 

* 

A12 = ( -0. 4480 )/( JGB) 

* 

A13 = ( 13. 62 4 ) /( JGB) 

* 

A21 = (2. 1131+.91786*( NGL/NGB ) )/( JDB ) 

* 

A22 = (-. 0998-0. 445*( NGL/NGB) )/( JDB) 

* 

A23 =( -0. 18851 )/( JDB) 

* 

A3 1 = 0. 0 

★ 

A3 2 = 0. 0 

* 

A33 = -1.0/TC 

★ 

Bll = 0. 0 

★ 

B12 * 0. 0 

★ 

B21 =. 035250/( JDB) 

★ 

B22 = -1. 0/( JDB) 

* 

B31 = 1.0/TC 

★ 

B32 = 0. 0 



DERIVATIVE 

NOSORT 

★ 

* COMPUTE INPUT TO THE NONLINEAR MODEL, MF( T) , WW(T). 

* RUN #9 

DMF =( MFDLY - MFI )/MFB 
DQD =( QDD - QDI)/QDB 

DNGDOT = All* DNG + A12*DNS + A13*DE + Bll *DMF 

DNSDOT = A2 1 * DNG + A22*DNS + A23*DE + B21*DMF + B22 *DQD 

DEDOT = A3 3 *DE + B31*DMF 



* DYNAMIC EQUATIONS FOR LINEAR MODEL. 
DNG=INTGRL( 0. 0, DNGDOT) 

DNS = I NTGRL( 0. 0,DNSDOT) 

DE =INTGRL( 0. 0, DEDOT) 

NGL = NGI + DNG* NGB 
NSL = NS I + DNS*NSB 

SORT 
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* THE STATEMENTS IN THE PREVIOUS (DERIVATIVE) SECTION YIELD VALUES 

* OF 'NG', AND 'NS' AS CALCULATED BY THE NONLINEAR AND STATE SPACE 

* MODELS. THE STATEMENTS BELOW COMPUTE THE VALUES OF 'NG', AND 'NS' 

* AS RECORDED FROM GAS TURBINE TEST DATA. 

* 

TERMINAL 
METHOD RK5 

CONTROL FINTIM=50. 0,DELT=0. 001 

PRINT 0. 5,NGD,NGL,NSD,NSL,MFD,QDD,DNSDOT,A21,DNG, . . . 

A22 , DNS , A23 , DE, B21 , DMF , B22 , DQD 

END 

STOP 



6 7 



APPENDIX E 



SOURCE DATA RUNS: 



Run 1 : Data run for dynamic Ns run at constant Ng. 

Ng — 20230 - 20240 rpm Ns — 505 - 2215 rpm 

Qd - 125 - 25 lbft Mf - 78.4 - 82.1 lbm/hr 

Run 3: Data run for dynamic Ng run at constant Ns. 

Ng — 20220 - 32180 rpm Ns -- 506 - 513 rpm 
Qd - 125 - 430 lbft Mf - 77.5 - 160.3 lbm/hr 

Run 5: Data run for dynamic Ng run at constant Ns 

Ng - 20220 - 32130 rpm Ns -- 1492 - 1510 rpm 

Qd - 71 - 144 lbft Mf- 83.8 - 165 lbm/hr 

Run 7 : Data run for dynamic Ng run at Ns (2252-2516 rpm). 

Ng - 20220 - 32170 rpm Ns - 2252 - 2516 rpm 

Qd - 25 - 240 lbft Mf - 78.8 - 170 lbm/hr 

Run 9: Data run for dynamic Ns run at constant Ng. 

Ng -- 32150 rpm Ns -- 510 - 2525 rpm 

Qd - 394 - 210 lbft Mf - 155.8 - 168.3 lbm/hr 
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